c  This program reads the unformatted head file from the model output
c  and makes a list of the computed heads versus measured heads.
c  The measured heads must be in a file that contains layer,row,column
c  water-level elevation, and any description desired.  The number of
c  known water levels is currently limited to 200. The format of
c  the file must be (i5,i5,i5,f10.2,a90)
c
c  modified and converted to DG by Lynette Brooks
c  September 11, 1991
c
      PARAMETER LENX=40000
      real X(20,200,200),value,head1(400),total,abstotal,errsumsq
      real head2(400),diff(400),tempx(200,200)
      integer arrow(400),arcol(400),arlay(400)
      real mean,absmean,resvar,dry
      integer count,maxcol,seqnum,rowcol
      CHARACTER FILENAME*96,ans*1,rest*90
      dimension text(4)
c  X array is the model head values for all model layers and the 
c    stress period, time step requested
c  value is the known vertical head diff between specified lays at a
c    specified row,col
c  diff is an array containing the difference between the modeled
c      value in the first and second specified layers
c  arrow,arcol,arlay are arrays of the layer, row, and column
c  maxcol is the maximum number of columns in the model
c  total is total of differences
c  abstotal is total of absolute difference
c  errsumsq is error sum of squares
c  mean is mean of differences
c  absmean is mean of absolute differences
c  resvar is residual variance
c  dry is the value printed if the cell dries up during model run
c  count counts the number of known water levels
c  rest is description of known water levels
c  text(4) is needed to read unformatted file correctly
c
C OPEN NODE INTEREST FILE
   10 write(6,*)' ENTER NAME OF FILE CONTAINING VERTICAL HEAD DIFFS: '
      read (5,'(A)',ERR=10) FILENAME
      OPEN (unit=1,FILE=FILENAME,STATUS='OLD',ERR=10)
C OPEN UNFORMATTED INPUT FILE
   20 write(6,*)' ENTER NAME OF UNFORMATTED MODEL HEAD OUTPUT FILE: '
      read (5,'(A)',ERR=20) FILENAME
      OPEN (unit=2,FILE=FILENAME,STATUS='OLD',FORM='UNFORMATTED',
     CERR=20)
      REWIND (unit=2)
C OPEN OUTPUT FILE
   30 write(6,*)' ENTER NAME OF OUTPUT FILE: '
      read (5,'(A)',ERR=30) FILENAME
      OPEN (UNIT=7,FILE=FILENAME,STATUS='UNKNOWN',ERR=30)
c
c
c
C READ FIRST ARRAY RECORD FOR PROBLEM DIMENSIONS
c
      read (2) KSTP,KPER,PERTIM,TOTIM,TEXT,NCOL,NROW,KLAY
C ALLOCATE SPACE IN THE X array and stop program if not large enough
      NODES=NCOL*NROW
      write (6,1010) nodes,LENX
      IF (LENX.lt.nodes) then
         write (6,*)' Model array is too large.  Program halts'
         STOP
      end if   
c
c  rewind file so it can be read again
c
      REWIND (unit=2)
c
c  initialize the counter for number of comparisons
c
c  count = 0
c
c
      istep = 0
      iper = 0
      write(6,*)' what time step,stress period,how many model layers?'
         read(5,*)istep,iper,ilay


c   
c  READ A LAYER OF DATA
c   read the first line of unformated head file
   50 read (2,END=70) KSTP,KPER,PERTIM,TOTIM,TEXT,NCOL,NROW,KLAY
      if (istep.ne.0) then
         if ((kstp.ne.istep).or.(kper.ne.iper)) then
             read (2,END=70) ((tempX(I,j),j=1,ncol),I=1,nrow)
c   reads the next line which are the heads 
c if the ts or str is not correct, read the next record
            go to 50
         else
            go to 60
         end if
      end if
c
c  if user hasn't specified only one set of data, prompt for each one
c
c  WRITE OUTPUT FILE FOR THIS LAYER if wanted
C  rewind the input file
      REWIND (unit=1)
c
c  READ IN THE COMPUTED HEADS FOR ALL LAYERS
c
 60   do 5,I = 1,ilay
        read(2,END=71) ((X(I,J,K),K=1,NCOL), J=1,NROW)
        read (2,END=71) KSTP,KPER,PERTIM,TOTIM,TEXT,NCOL,NROW,KLAY
 5    CONTINUE
 71   CONTINUE
      WRITE(7,*)' UPPER LAY IS SUBTRACTED FROM LWR LAY,
     1 + IS UP VERT GRAD'
 63   write(7,1030)kper,kstp,klay
      write(7,*)' ROW  COL  KNOW-DIFF  COMP-DIFF CALC-DIFF  HEAD1
     1    HEAD2        DESCRIPTION'
c
c  read known head values one line at a time, calculate difference
c  and write data
c
   62 read (1,'(4I5,f10.2,a90)',END=70) LAY1,LAY2,IROW,ICOL,VALUE,rest
c   read the input file, only lines that correspond to the KLAY
      count = count + 1
      if (x(lay1,irow,icol).gt.1000000.or.x(lay1,irow,icol).lt.0.0001)
     C then
         head1(count)= dry
      else
         head1(count) = x(LAY1,irow,icol)
         HEAD2(COUNT) = x(LAY2,irow,icol)
         TEMP = HEAD2(COUNT) - HEAD1(COUNT)
         DIFF(COUNT)  = TEMP - VALUE
      end if
      arrow(count)=irow
      arcol(count)=icol
      arlay(count)=lay

c
      WRITE (7,'(2I5,5f10.2,a90)') IROW,ICOL,VALUE,
     1        temp,diff(count),HEAD1(COUNT),HEAD2(COUNT),rest
      GO TO 62
c
c
c  all data read, do statistics
c  initialize totals
   70 total = 0.
c      abstotal = 0.
c      errsumsq = 0.
c      write(6,*)' Number of Comparisons =         ',count
cc
cc  calculate totals
cc
c      do 80 i =1,count
c         total = total + diff(i)
c         abstotal = abstotal + abs(diff(i))
c         errsumsq = errsumsq + (diff(i)**2)
c 80   continue
cc
cc  calculate statistics and write results
cc
c      mean = total/count
c      absmean = abstotal/count
c      resvar = errsumsq/(count-2)
cc
c      write(7,*)' Number of Comparisons =         ',count
c      write(7,*)' Total of Differences =          ',total
c      write(7,*)' Total of Absolute Differences = ',abstotal
c      write(7,*)' Mean of Differences =           ',mean
c      write(7,*)' Mean of Absolute Differences =  ',absmean
c      write(7,*)' Error Sum of Squares =          ',errsumsq
c      write(7,*)' Residual Variance =             ',resvar
cc
cc  prompt for and open ascii file that can be read by info
cc
c      write(6,*)'do you want an ascii file for info?'
c      read(5,'(a)')ans
c      if (ans.eq.'y'.or.ans.eq.'Y') then
c
c         write(6,*)'enter name of ascii info file'
c         read(5,*)filename
c         open (8,file=filename,status='unknown')
c
c  prompt for number of columns and do calculations and printout
c
c         write(6,*)'enter number of columns in model'
c         read(5,*)maxcol
c         do 90 i=1,count
c            seqnum = (arrow(i)-1)*maxcol + arcol(i)
c            rowcol = 100*arrow(i)+arcol(i)
c            idiff  = diff(i)
c            write(8,1040)seqnum,rowcol,arlay(i),diff(i),idiff
c 90      continue
c      end if
 1010 FORMAT (I10,' ELEMENTS OF THE X ARRAY USED OUT OF',I10,
     -' POSSIBLE.')
 1020 format(' do you want data for time step',i5,' period',i5,
     1           ' and layer',i5,'(y/n)?')
 1030 format(' Data for Stress Period',i5,' Time Step',i5,
     1            ' Layer',i5)
 1040 format (i5,',',i5,',',i5,',',f10.2,','i10)
cc
cc  close files and stop
cc
      CLOSE (unit=1)
      CLOSE (unit=2)
      ENDFILE (UNIT=7)
      CLOSE (UNIT=7)
      stop
      END
